%!Tex Program = xelatex
\documentclass[a4paper]{ctexart}
\usepackage{xltxtra}
\usepackage{listings}
\usepackage{float}
\usepackage{graphicx}
\usepackage{geometry}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{enumitem}
\geometry{a4paper,scale=0.75}
\title{Step-3 Tutorial Introduction}
\author{周川迪}
\date{\today}
\begin{document}
\maketitle
\subsection*{有限元方法的基本设置}

这是我们实际使用有限元方法计算某些东西的第一个例子。我们将求解一个简化版本的泊松方程，边界值为零，右手边为非零：
 \begin{align*} -\Delta u &= f \qquad\qquad & \text{in}\ \Omega, \\ u &= 0 \qquad\qquad & \text{on}\ \partial\Omega. \end{align*} 

我们将在正方形 $\Omega=[0,1]^2$,上求解这个方程，你已经学会了如何在步骤1和步骤2中生成网格。在这个程序中，我们也只考虑特殊情况$f(\mathbf x)=1$ ，并在下一个教程程序（步骤4）中讨论如何实现更一般的情况。

如果你已经了解了有限元方法的基础知识，你会记得我们需要采取哪些步骤来用有限维逼近来近似求解$u$。具体来说，我们首先需要导出上述方程的弱形式，这是通过从左边用一个测试函数$\varphi$乘以方程，并在域$Ω$上积分得到的（我们将在下面回到从左边而不是从右边乘法的原因）：
\begin{align*} -\int_\Omega \varphi \Delta u = \int_\Omega \varphi f. \end{align*} 
这可以通过分部积分得出：
\begin{align*} \int_\Omega \nabla\varphi \cdot \nabla u - \int_{\partial\Omega} \varphi \mathbf{n}\cdot \nabla u = \int_\Omega \varphi f. \end{align*} 
测试函数$\varphi$必须满足与边界条件的相同类型（在数学术语中：它需要来自我们寻找解的集合的切空间），因此在边界上$\varphi=0$，因此我们要找的弱形式变为 
\begin{align*} (\nabla\varphi, \nabla u) = (\varphi, f), \end{align*} 

这里我们使用通常的表示法 $(a,b)=\int_\Omega a\; b$。然后问题要求找到一个函数$u$，对于所有适当空间中的测试函数$\varphi$，都满足这个陈述。

当然，在一般情况下，我们无法在计算机上找到这样一个函数，
而是寻求一个近似的解$u_h(\mathbf x)=\sum_j U_j \varphi_j(\mathbf x)$，其中$U_j$是我们需要确定的未知展开系数(问题的"自由度")，$\varphi_i(x)$是我们将使用的有限元形状函数。为了定义这些形状函数，我们需要以下内容：
\begin{itemize}[label=$\bullet$]
    \item 在其上定义形状函数的网格。你已经看到了如何在步骤1和步骤2中生成和处理描述网格的对象。
    \item 描述我们想要在参考单元（在deal.II中始终是单位区间$[0,1]$、单位正方形$[0,1]^2$或单位立方体$[0,1]^3$，取决于你所工作的空间维数）上使用的形状函数的有限元。在步骤2中，我们已经使用了类型为FE\_Q<2>的对象，它表示以支持点上的插值定义形状函数的常规拉格朗日元素。最简单的是FE\_Q<2>(1)，它使用多项式次数为1。在二维中，它们通常被称为双线性，因为它们在参考单元的两个坐标中的每一个上都是线性的。（在一维中，它们是线性的，在三维中是三线性的。然而，在deal.II的文档中，我们经常不进行这种区分，只是简单地称这些函数为"线性"。）
    \item DoFHandler对象，在网格上枚举所有自由度，以参考单元描述有限元对象提供的基础。你还在步骤2中看到了如何做到这一点。
    \item 一个映射，告诉如何从参考单元上有限元类定义的形状函数获得实际单元上的形状函数。默认情况下，除非你明确说明，deal.II将为此使用一个（双线性、三线性）映射，所以在大多数情况下，你不必担心这一步。
\end{itemize}

通过这些步骤，我们现在有一组函数$\varphi_i$，我们可以定义离散问题的弱形式：找到一个函数$u_h$,，即找到上面提到的展开系数 $U_j$ ，使得
\begin{align*} (\nabla\varphi_i, \nabla u_h) = (\varphi_i, f), \qquad\qquad i=0\ldots N-1. \end{align*} 

请注意，我们在这里遵循的约定是一切都从零开始计数，与C和C++中常见的约定相同。如果你插入表示 $u_h(\mathbf x)=\sum_j U_j \varphi_j(\mathbf x)$ ，会观察到
\begin{align*} (\nabla\varphi_i, \nabla u_h) &= \left(\nabla\varphi_i, \nabla \Bigl[\sum_j U_j \varphi_j\Bigr]\right) \\ &= \sum_j \left(\nabla\varphi_i, \nabla \left[U_j \varphi_j\right]\right) \\ &= \sum_j \left(\nabla\varphi_i, \nabla \varphi_j \right) U_j. \end{align*}
有了这个，问题可以重新表述为一个线性系统：找到一个向量$U$，使得
\begin{align*} A U = F, \end{align*} 
其中矩阵$A$和右手边$F$定义为
\begin{align*} A_{ij} &= (\nabla\varphi_i, \nabla \varphi_j), \\ F_i &= (\varphi_i, f). \end{align*} 

\subsection*{我们应该从左边还是从右边乘以一个测试函数？}

在我们继续描述如何计算这些量之前，注意如果我们将原始方程从右边而不是从左边乘以一个测试函数，那么我们将获得一个形式为
\begin{align*} U^T A = F^T \end{align*} 
的线性系统，其中$F^T$是一个行向量。通过转置这个系统，显然等价于求解
\begin{align*} A^T U = F \end{align*}
这在这里与上面的情况相同，因为$A=A^T$。但一般来说不是这样的，在避免任何混淆的目的，经验表明，习惯上从左边而不是从右边乘以方程（如数学文献中经常这样做）能避免一类常见错误，因为矩阵在理论和实现之间自动正确，无需转置。参见本教程的步骤9，这是我们第一个例子，在这个例子中，我们有一个非对称的双线性型，从右边乘以与从左边乘以有所不同。

\subsection*{计算矩阵和右手边向量}

在计算矩阵和右手向量之前，我们需要以下几个对象和计算方法：

\begin{itemize}[label=$\bullet$]
    \item 对于矩阵$A$、向量$U$和向量$F$，它们的类型是稀疏矩阵和向量。在下面的程序中，我们将看到哪些类用于解线性系统。
    \item 我们需要一种计算积分的方法。在有限元方法中，最常用的方法是使用数值积分，即将积分替换为对每个单元上一组点的加权求和。首先，我们将整个域$\Omega$上的积分拆分为各个单元上的积分：
    \begin{align*}
    A_{ij} = (\nabla\varphi_i, \nabla \varphi_j) = \sum_{K \in {\mathbb T}} \int_K \nabla\varphi_i \cdot \nabla \varphi_j, \
    F_i = (\varphi_i, f) = \sum_{K \in {\mathbb T}} \int_K \varphi_i f,
    \end{align*}
    然后通过数值积分来近似每个单元的贡献：
    \begin{align*}
    A^K_{ij} = \int_K \nabla\varphi_i \cdot \nabla \varphi_j \approx \sum_q \nabla\varphi_i(\mathbf x^K_q) \cdot \nabla \varphi_j(\mathbf x^K_q) w_q^K, \
    F^K_i = \int_K \varphi_i f \approx \sum_q \varphi_i(\mathbf x^K_q) f(\mathbf x^K_q) w^K_q,
    \end{align*}
    其中 $\mathbf x^K_q$ 是单元$K$上的第$q$个数值积分点，$w^K_q$是第$q$个数值积分权重。在进行这个过程时，需要完成不同的任务，下面我们将逐个讨论。
    \item 首先，我们需要一种描述数值积分点位置 $\mathbf x_q^K$及其权重 $w^K_q$的方法。通常，这些点和权重是从参考单元映射到实际单元上，也就是使用 MappingQ1 类隐式地完成映射，或者如果你明确指定了其他类，则通过继承自 Mapping 的其他类来完成。描述参考单元上的积分点和权重的对象是从 Quadrature 基类派生的。通常，我们选择一个数值积分公式（即一组点和权重），使得数值积分恰好等于矩阵中的积分；这可以通过使用高斯数值积分公式来实现，其中 QGauss 类提供了高斯数值积分的实现。
    \item 接下来，我们需要一种帮助我们在单元$K$上计算 $\varphi_i(\mathbf x^K_q)$的方法。这就是 FEValues 类做的事：它接受一个有限元对象来描述参考单元上的$\varphi$，一个数值积分对象来描述数值积分点和权重，以及一个映射对象（或者隐式地使用 MappingQ1 类），并在实际单元$K$上的数值积分点上提供形状函数和导数的值以及积分所需的其他信息。
\end{itemize}

FEValues 类真正是装配过程中的核心类。你可以将其视为以下方式：FiniteElement 和其派生类描述形状函数，即无限维对象：函数在每个点都有值。出于理论原因，我们需要这样做，因为我们希望用函数的积分来进行分析。然而，对于计算机来说，这是一个非常困难的概念，因为它们通常只能处理有限量的信息，因此我们通过数值积分点的求和来代替积分，这些积分点是通过映射（Mapping 对象）使用参考单元（Quadrature 对象）上定义的点映射到实际单元上的。本质上，我们将问题简化为只需要有限量的信息，即形状函数的值和导数、数值积分权重、法向量等，仅在有限一组点上。FEValues 类是将这三个组件集成在一起并提供特定单元$K$上的有限信息集的类。在下面组装线性系统时，你将看到它的作用。

值得注意的是，如果你在应用程序中自己创建这三个对象，并自己处理信息，也可以实现同样的功能。然而，这既不会更简单（FEValues 类提供了你实际需要的确切信息），也不会更快：FEValues 类经过高度优化，只计算你在每个单元上实际需要的特定信息；如果上一个单元的信息可以重复使用，则会重复使用，该类中有很多代码确保在有利时进行缓存。

最后一点是，在获得线性系统之后，使用迭代求解器求解，然后进行后处理：我们使用 DataOut 类创建输出文件，然后可以使用常见的可视化程序对其进行可视化。

\subsection*{关于应用}
尽管这是使用有限元法求解的最简单方程，但该程序展示了大多数有限元程序的基本结构，并且也作为后续几个程序的模板。具体来说，该程序的主函数如下所示：
\begin{lstlisting}[language=C++]
class Step3
{
  public:
    Step3 ();
    void run ();
 
  private:
    void make_grid ();
    void setup_system ();
    void assemble_system ();
    void solve ();
    void output_results () const;
 
    Triangulation<2>     triangulation;
    FE_Q<2>              fe;
    DoFHandler<2>        dof_handler;
 
    SparsityPattern      sparsity_pattern;
    SparseMatrix<double> system_matrix;
    Vector<double>       solution;
    Vector<double>       system_rhs;
};
\end{lstlisting}

我们从成员变量开始：这些成员变量遵循我们在上面的要点中概述的构建模块，即我们需要一个Triangulation和一个DoFHandler对象，以及描述我们要使用的形状函数类型的有限元对象。第二组对象与线性代数有关：系统矩阵和右手边以及解向量，以及描述矩阵稀疏模式的对象。这就是该类需要的所有内容（以及任何一个解决器用于静态PDE的必要要素），并且需要在整个程序中保持有效。与此相反，我们需要用于组装的FEValues对象仅在组装过程中需要，因此我们在执行组装的函数中将其创建为局部对象，并在其结束时销毁。

再看一下成员函数。这些函数已经形成了几乎所有后续教程程序将使用的常见结构：
\begin{itemize}[label=$\bullet$]
    \item make\_grid(): 这可以称为预处理函数。顾名思义，它设置了存储三角网格的对象。在后面的示例中，它还可以处理边界条件、几何形状等。
    \item setup\_system(): 这是设置解决问题所需的所有其他数据结构的函数。特别是，它将初始化DoFHandler对象并正确调整与线性代数相关的各种对象的大小。将此函数与上面的预处理函数分开是因为，在时间相关的程序中，每当网格自适应细化（我们将在步骤6中学习如何做到这一点）时，它可能至少被调用几次。另一方面，在上面的预处理函数中设置网格仅在程序开始时执行一次，因此将其分离成自己的函数。
    \item assemble\_system(): 这是计算矩阵和右手边的内容的地方，正如上面的介绍中详细讨论的那样。由于对该线性系统进行某些操作在概念上与计算其条目非常不同，因此我们将其与以下函数分开。
    \item solve(): 这是计算线性系统$AU=F$的解$U$的函数。在当前程序中，由于矩阵非常简单，因此这是一个简单的任务，但是当问题不再那么简单时（参见步骤20、步骤22或学习了有关该库的更多内容后，在步骤31中），计算解的函数将成为程序大小的重要部分。
    
    \item output\_results(): 最后，当您计算出解时，可能希望对其进行某些操作。例如，您可能希望以可视化格式输出它，或计算您感兴趣的量：例如热交换器中的热通量、翼的空气摩擦系数、最大桥载荷，或者仅在某一点的数值解的值。因此，此函数是进行后处理的地方。
\end{itemize}
所有这些都由单个公共函数（除构造函数外）run()函数组成。从创建此类型对象的位置调用该函数，并且该函数按正确的顺序调用所有其他函数。将这个操作封装到run()函数中，而不是从main()函数中调用所有其他函数，可以确保您可以更改该类内部实现的关注点分离方式。例如，如果某个函数变得太大，您可以将其拆分成两部分，作为结果唯一需要关注更改的位置就是该类内部，而不是任何其他地方。

在deal.II中，通过在types命名空间中使用别名定义了一些以 types::global\_dof\_index 的形式出现的整数类型，用于表示自由度的全局索引，即在DoFHandler对象中的特定自由度的索引（而不是在特定单元内的特定自由度的索引）。对于当前程序（以及几乎所有教程程序），全局未知数通常有几千到几百万个（对于$Q_1$元素，在二维中每个单元上有4个本地自由度，在三维中有8个）。因此，unsigned int 数据类型足以存储全局自由度索引，因为它允许存储从0到略多于40亿的数字（在大多数系统上，整数是32位）。

那么，为什么不直接使用 unsigned int 呢？deal.II 在7.3版本之前确实是这样做的。然而，deal.II 支持非常大的计算（通过 step-40 中讨论的框架），当跨越几千个处理器时，可能会有超过40亿个未知数的情况。因此，在某些情况下，unsigned int 不足够大，我们需要一个64位的无符号整数类型。为了实现这一点，我们引入了 types::global\_dof\_index，默认情况下它被定义为 unsigned int，但可以通过在配置期间传递特定标志来将其定义为 unsigned long long int（请参阅ReadMe文件）。

在库和构建在其上的代码中，如果看到使用 types::global\_dof\_index 数据类型的地方，您立即就知道所引用的数量实际上是全局自由度索引。如果我们只是使用 unsigned int（它也可能是一个本地索引、边界指示器、材料ID等），则不会出现这样的含义。立即知道变量所引用的内容也有助于避免错误：如果您看到将 types::global\_dof\_index 类型的对象赋给 types::subdomain\_id 类型的变量，尽管它们都由无符号整数表示且编译器不会报错，但很明显存在错误。

从更实际的角度来说，存在这种类型的存在意味着在组装过程中，我们会创建一个 $4\times4$ 的矩阵（在二维中使用$Q_1$元素），这个矩阵包含当前单元的贡献，并且我们需要将这个矩阵的元素添加到全局（系统）矩阵的相应元素中。为此，我们需要获得局部于当前单元的全局自由度的全局索引，我们将始终使用以下代码片段：
\begin{lstlisting}[language=C++]
cell->get_dof_indices (local_dof_indices);
\end{lstlisting}
其中 local\_dof\_indices 被声明为
\begin{lstlisting}[language=C++]
std::vectortypes::global_dof_index //接下面一行
local_dof_indices(fe.n_dofs_per_cell());
\end{lstlisting}
这个变量的名称可能有点误导性，它代表着“当前单元上局部定义的自由度的全局索引”，但是在整个库中，持有这些信息的变量普遍以这种方式命名。

\subsection*{关于step-3计算的问题的描述}
deal.II是一个用于求解偏微分方程的开源有限元库。在deal.II的examples中，step-3是一个用于求解二维Poisson方程的示例程序。

step-3程序使用有限元方法来解决Poisson方程的边界值问题。Poisson方程是一个常见的椭圆型偏微分方程，常用于描述许多物理和工程问题。它的一般形式可以表示为：
\begin{align*}
    \nabla^2 u=f
\end{align*}
其中，$u$是未知函数，$\nabla^2$是拉普拉斯算子，$f$是给定的函数。

具体参数如下：
\begin{itemize}[label=$\bullet$]
    \item 网格：step-3使用一个简单的二维正方形网格，可以通过设置网格的大小来调整分辨率。
    \item 有限元空间：程序使用连续Galerkin有限元空间来近似解。你可以通过设置多项式的次数来调整有限元的精度。
    \item 有限元方法：程序使用标准的拉格朗日有限元方法，使用三角形或四边形单元作为离散化的基本单元。
    \item 边界条件：程序指定了Dirichlet边界条件，即在边界上给定了特定的函数值。
\end{itemize}

运行step-3程序后，它会生成一个名为solution.gnuplot的数据文件，其中包含了近似解的节点坐标和数值。使用gnuplot可将这些数据绘制成图像如下：

\begin{figure}[h!]
    \centering
    \includegraphics{solution.png} 
\end{figure}

\end{document}